act_site = [1 2 4 7 8];
pas_site = [3 5 6 9];
con_site = 10;

for site = 1:9
    for layer = 1:4
        for time = 1:4
            cur_fb_act = proj_meta(site).rd(layer,time).act(:,1:7500);
            cur_pb_act = proj_meta(site).rd(layer,time).act(:,7501:15000);

            % generate binary activity during feedback.
            cur_fb_bi = -proj_meta(site).rd(layer,time).velM_smoothed(1:7500);
            cur_fb_bi(cur_fb_bi > 0.001) = 1;
            cur_fb_bi(cur_fb_bi <= 0.001) = 0;
            
            % generate binary activity during playback.
            cur_pb_bi = -proj_meta(site).rd(layer,time).velM_smoothed(7501:15000);
            cur_pb_bi(cur_pb_bi > 0.001) = 1;
            cur_pb_bi(cur_pb_bi <= 0.001) = 0;
            
            % generate the frame id of ps onset.
            ps_ons = find(diff(proj_meta(site).rd(layer,time).ps_id(1:7500))>1)+1;
            % remove early & late ps.
            ps_ons(ps_ons > 7425) = [];
            ps_ons(ps_ons < 6) = [];
            
            % generate the frame id of playback halt onset.
            pb_ons = find(diff(proj_meta(site).rd(layer,time).ps_id(7501:15000))>1)+1;
            % remove early & late pb.
            pb_ons(pb_ons > 7425) = [];
            pb_ons(pb_ons < 8) = [];

            rps_nbr = 0; % number of real ps
            rps_id = [];
            for ps_id = 1:length(ps_ons)
                if sum(cur_fb_bi(ps_ons(ps_id)-5:ps_ons(ps_id)+15)) == 21
                    rps_nbr = rps_nbr+1;
                    rps_id(rps_nbr) = ps_id;
                    all_rps_nbr(site,time) = rps_nbr;
                end
            end
            
            rpb_nbr = 0; % number of real pb halt
            rpb_id = [];
            for pb_id = rps_id % rpb is selected based on rps
                if sum(cur_pb_bi(pb_ons(pb_id)-7:pb_ons(pb_id)+22)) == 0
                    rpb_nbr = rpb_nbr+1;
                    rpb_id(rpb_nbr) = pb_id;
                    all_rpb_nbr(site,time) = rpb_nbr;
                end
            end

            if rps_nbr >= 3
                rps_id = rps_id(1:3);
            else
                rps_id = [];
            end
            
            if rpb_nbr >= 3
                rpb_id = rpb_id(1:3);
            else
                rpb_id = [];
            end
            
            for cell_id = act_cell{site,layer}
                rps_resp_cnt = 0;
                prps_resp = [];
                fire_fb = 0;
                siln_fb = 0;
                for rps_act_nbr = 1:length(rps_id)
                    rps_frame_id = ps_ons(rps_id(rps_act_nbr));
                    prps_win = rps_frame_id+8:rps_frame_id+22;
                    rps_resp_cnt = rps_resp_cnt+1;
                    if max(cur_fb_act(cell_id,prps_win)) > crit{site,layer,cell_id}
                        prps_resp(rps_resp_cnt) = max(cur_fb_act(cell_id,prps_win))-1;
                        fire_fb = fire_fb+1;
                    else
                        prps_resp(rps_resp_cnt) = 0;
                        siln_fb = siln_fb+1;
                    end
                end
                prps_resp_by_cell{site,layer,time,cell_id} = prps_resp;
                fire_fb_cnt{site,layer,time,cell_id} = fire_fb;
                siln_fb_cnt{site,layer,time,cell_id} = siln_fb;
            end
            
            for cell_id = act_cell{site,layer}
                rpb_resp_cnt = 0;
                prpb_resp = [];
                fire_pb = 0;
                siln_pb = 0;
                for rpb_act_nbr = 1:length(rpb_id)
                    rps_frame_id = ps_ons(rpb_id(rpb_act_nbr));
                    prpb_win = rps_frame_id+8:rps_frame_id+22;
                    rpb_resp_cnt = rpb_resp_cnt+1;
                    if max(cur_pb_act(cell_id,prpb_win)) > crit{site,layer,cell_id}
                        prpb_resp(rpb_resp_cnt) = max(cur_pb_act(cell_id,prpb_win))-1;
                        fire_pb = fire_pb+1;
                    else
                        prpb_resp(rpb_resp_cnt) = 0;
                        siln_pb = siln_pb+1;
                    end
                end
                prpb_resp_by_cell{site,layer,time,cell_id} = prpb_resp;
                fire_pb_cnt{site,layer,time,cell_id} = fire_pb;
                siln_pb_cnt{site,layer,time,cell_id} = siln_pb;
            end
        end
        cell_cnt{site,layer} = size(proj_meta(site).rd(layer,1).act,1);
    end
end

for site = 1:9
    for layer = 1:4
        if act_cell_cnt{site,layer} > 0
            for cell_id = act_cell{site,layer}
                cur_cell_fb_act = [prps_resp_by_cell{site,layer,:,cell_id}];
                mprps{site,layer,cell_id} = mean(cur_cell_fb_act);
                
                fire_fb_cnt_s = sum([fire_fb_cnt{site,layer,:,cell_id}]);
                siln_fb_cnt_s = sum([siln_fb_cnt{site,layer,:,cell_id}]);
                prps_prob{site,layer,cell_id} = fire_fb_cnt_s/(fire_fb_cnt_s+siln_fb_cnt_s);
            end
        end
    end
end

for site = 1:9
    for layer = 1:4
        if act_cell_cnt{site,layer} > 0
            for cell_id = act_cell{site,layer}
                cur_cell_pb_act = [prpb_resp_by_cell{site,layer,:,cell_id}];
                mprpb{site,layer,cell_id} = mean(cur_cell_pb_act);
                
                fire_pb_cnt_s = sum([fire_pb_cnt{site,layer,:,cell_id}]);
                siln_pb_cnt_s = sum([siln_pb_cnt{site,layer,:,cell_id}]);
                prpb_prob{site,layer,cell_id} = fire_pb_cnt_s/(fire_pb_cnt_s+siln_pb_cnt_s);
            end
        end
    end
end